function NC = water_extracted(f,f_watered,watermark,varargin)
% 水印提取 + 效果检测函数
% 
% Input:
%   f:原始未加水印图像
%   f_watered:已添加水印图像
%   watermark:水印图像
% 
% Output:
%   NC:归一化交叉相关（Normalized Cross-Correlation，NC）系数
%   
%   ......
% 
    if ~isempty(varargin)
        alpha = varargin{1};
    else
        alpha = 0.1;
    end


    % f_watered = SURF_Correct(f,f_watered);

    % 含水印图像四尺度变换
    y = nsctdec(f,[2,4],'pkva','9-7');
    y_watered = nsctdec(f_watered,[2,4],'pkva','9-7');
    
    % svd分解提取水印信息
    [U_original1,S_original1,V_original1] = svd(y{1});
    [U_original2,S_original2,V_original2] = svd(y{3}{5});
    [U_watered1,S_watered1,V_watered1] = svd(y_watered{1});
    [U_watered2,S_watered2,V_watered2] = svd(y_watered{3}{5});
    [U,S,V] = svd(watermark);
    % 
    water_extracted_S = zeros(64,64);

    % alpha = 0.1; % 嵌入强度设置为0.2
    water_extracted_S1 = (S_watered1(1:32,1:32) - S_original1(1:32,1:32)) / (alpha);
    % water_extracted_S1 = (S_watered1(1:32,1:32)) / (0.05);
    water_extracted_S2 = (S_watered2(1:32,1:32) - S_original2(1:32,1:32)) / (0.2);

    water_extracted_S(1:32,1:32) = water_extracted_S1;
    water_extracted_S(33:64,33:64) = water_extracted_S2;

    % water_extracted = zeros(64,64);
    % for i = 1:8
    %     for j = 1:8
    %         [U_y,S_y,V_y] = svd(y{1}(i*8-7:i*8,j*8-7:j*8));
    %         [U_w,S_w,V_w] = svd(y_watered{1}(i*8-7:i*8,j*8-7:j*8));
    %         a((i-1)*8+j) = (S_w(1,1)-S_y(1,1))/0.2;
    %     end
    % end
    % water_extracted_S = diag(a);

    % % svd分解提取水印信息
    % [U_original,S_original,V_original] = svd(y{3}{5});
    % [U_watered,S_watered,V_watered] = svd(y_watered{1});
    % [U,S,V] = svd(watermark);
    % 
    % alpha = 0.1; % 嵌入强度设置为0.2
    % water_extracted_S = (S_watered(1:128,1:128) - S_original(1:128,1:128)) / (alpha);


    % 水印恢复，逆svd
    water_extracted = U * water_extracted_S * V';

    figure;imshow(water_extracted,[]);
    title("恢复的水印");

    % NC系数计算
    mu_water_extracted = mean(water_extracted(:));
    mu_watermark = mean(watermark(:));

    sigma_water_extracted = std(water_extracted(:));
    sigma_watermark = std(watermark(:));

    water_extracted_n = (water_extracted - mu_water_extracted) / sigma_water_extracted;
    watermark_n = (watermark - mu_watermark) / sigma_watermark;

    NC = sum(water_extracted_n(:) .* watermark_n(:)) / numel(water_extracted);
end